﻿#include <xy/math/poly.h>
#include <xy/solver/diff.h>

namespace xy
{

namespace solver
{

using namespace xy::math;

void ode::solve(funcxxt f, time_grid &G, std::shared_ptr<multi_step> Z)
{
    // 指针位置
    gtimer U = G.begin();

    std::size_t s = Z->series;
    std::size_t N = G.size();

    // Classical RK4 计算初步
    auto RK = std::make_shared<classical_RK4>(Z->stepTime);
    for (std::size_t i = 0; i < s - 1; i++)
        U = RK->shift(f, U);

    // 同步时间
    Z->currentTime = RK->currentTime;

    // 循环计算每一步，N 个节点，移动 N-1 次
    for (std::size_t i = 0; i < N - s; i++)
        U = Z->shift(f, U);
}

void ode::solve(funcxxt f, time_grid &G, std::shared_ptr<one_step> Z, bool extra)
{
    // 指针位置
    gtimer U = G.begin();

    // 循环计算每一步，N 个节点，移动 N-1 次
    std::size_t N = G.size();
    for (std::size_t i = 1; i < N; i++)
    {
        auto V = Z->shift(f, U);
        U = extra ? Z->extrapolate(f, U) : V;
    }
}

void pde::solve(const matxd &A, func1x g, time_grid &G, std::shared_ptr<ESDIRK6> Z)
{
    // 指针位置
    gtimer U = G.begin();

    // 系数矩阵 B = (I - k a_ii A) 由于对角元相同，因此不重复计算
    auto B = A * A(0, 0) * Z->stepTime;
    for (std::size_t i = 0; i < A.rows(); i++)
        B(i, i) = 1 - B(i, i);

    // 线性求解器
    linear_ptlu solver(std::move(B));
    solver.perform();

    // 循环计算每一步，N 个节点，移动 N-1 次
    std::size_t N = G.size();
    for (std::size_t i = 1; i < N; i++)
        U = Z->shift(A, g, solver, U);
}

} // namespace solver

} // namespace xy
